#!/bin/tcsh -fe
#############################################################################
#                                                                           #
# Title: CTFFIND4.1                                                         #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 09/39/2014                                             #
# Last Modification: 03/08/2015	                                            #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 35
#
# MANUAL: This script measures the CTF of the recorded image, using CTFFIND4
#
# DISPLAY: imagenumber
# DISPLAY: comment
# DISPLAY: defocus
# DISPLAY: phacon
# DISPLAY: CS
# DISPLAY: KV
# DISPLAY: DEFOCUS_TLTAXIS
# DISPLAY: DEFOCUS_TLTANG
# DISPLAY: sample_pixel
# DISPLAY: phacon
# DISPLAY: CTFFIND4_defocus
# DISPLAY: defocus
# DISPLAY: df_start
# DISPLAY: df_end
# DISPLAY: df_step
# DISPLAY: CTFFIND4_defocus_res_min
# DISPLAY: CTFFIND4_defocus_res_max
# DISPLAY: defocus_CCvalue
# DISPLAY: defocus_RESMAX
# DISPLAY: defocus_phase_shift_doit
# DISPLAY: defocus_phase_shift
# DISPLAY: defocus_astig
# DISPLAY: defocus_defocus
# DISPLAY: CTFFIND4_phase_shift
# DISPLAY: CTFFIND4_defocus_phase_shift_L
# DISPLAY: CTFFIND4_defocus_phase_shift_H
# DISPLAY: CTFFIND4_defocus_phase_shift_step
# DISPLAY: CTFFIND4_inputfile
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
set app_2dx_mrc_converter = ""
#
set tempkeep = ""
set imagename = ""
set nonmaskimagename = ""
set imagenumber = ""
set imagesidelength = ""
set magnification = ""
set stepdigitizer = ""
set Calc_from_sample_pixel = ""
set sample_pixel = ""
set phacon = ""
set RESMIN = ""
set RESMAX = ""
set RADLIM = ""
set CS = ""
set KV = ""
set defocus = ""
set CTFFIND4_defocus = ""
set phacon = ""
set movie_stackname = ""
set CTFFIND4_RESMAX = ""
set df_start = ""
set df_end = ""
set df_step = ""
set CTFFIND4_defocus_res_min = ""
set CTFFIND4_defocus_res_max = ""
set defocus_RESMAX = ""
set defocus_CCvalue = ""
set defocus_phase_shift_doit = ""
set import_original_time = ""
set defocus_phase_shift = ""
set defocus_phase_shift_doit = ""
set defocus_astig = ""
set defocus_defocus = ""
set CTFFIND4_phase_shift = ""
set CTFFIND4_defocus_phase_shift_L = ""
set CTFFIND4_defocus_phase_shift_H = ""
set CTFFIND4_defocus_phase_shift_step = ""
set CTFFIND4_inputfile = ""
set raw_gaincorrectedstack = ""
#
#$end_vars
#
set scriptname = process_ctffind
\rm -f LOGS/${scriptname}.results
#
source ${proc_2dx}/initialize
#
echo "<<@progress: 5>>"
#
if ( ${defocus_phase_shift_doit} == "-" ) then
  set defocus_phase_shift_doit = ${Default_phase_shift_doit}
  echo "set defocus_phase_shift_doit = ${defocus_phase_shift_doit}" >> LOGS/${scriptname}.results
endif
#
if(${import_original_time} == "-" || ${import_original_time} == "") then
  @ status_date = `date +%s` * 1000
  set date_text = "Processed at "`echo ${status_date} | awk ' { s = $1 / 1000 } END { print s } ' | ${app_gawk} '{print strftime("%c", $0)}' `
else
  set status_date = ${import_original_time}
  set date_text = "Recorded at "`echo ${status_date} | awk ' { s = $1 / 1000 } END { print s } ' | ${app_gawk} '{print strftime("%c", $0)}' `
endif
#
set defocus_wrong = `echo ${CTFFIND4_defocus_res_min} 30.0 | awk '{ if ( $1 > $2 ) { s = 1 } else { s = 0 } } END { print s } '`
if ( ${defocus_wrong} == "1" ) then
  set CTFFIND4_defocus_res_min = "30.0"
  ${proc_2dx}/linblock "WARNING: correcting CTFFIND4_defocus_res_min to ${CTFFIND4_defocus_res_min}" 
  echo "set CTFFIND4_defocus_res_min = ${CTFFIND4_defocus_res_min}" >> LOGS/${scriptname}.results
endif
set ampcon = ` echo "scale=3; sqrt( 1 - ${phacon} * ${phacon} )" | bc `
#
if ( ${CTFFIND4_inputfile} == 0 && -e ${raw_gaincorrectedstack}.mrcs ) then
  set input_image = ${raw_gaincorrectedstack}
  set is_movie = "yes"
  ${proc_2dx}/linblock "Using raw gain-corrected stack ${input_image}.mrcs as input file for CTF measurement."
# else
#   ${proc_2dx}/protest "ERROR: Gain-corrected stack ${raw_gaincorrectedstack}.mrcs does not exist."
endif
if ( ${CTFFIND4_inputfile} == 1 && -e  ${movie_stackname}.mrcs ) then
  set input_image = ${movie_stackname}
  set is_movie = "yes"
  ${proc_2dx}/linblock "Using drift-corrected stack ${input_image}.mrcs as input file for CTF measurement."
# else
#   ${proc_2dx}/protest "ERROR: Drift-corrected stack ${movie_stackname}.mrcs does not exist."
endif
if ( ${CTFFIND4_inputfile} == 2 && -e  ${movie_stackname}_Sum.mrc ) then
  set input_image = ${movie_stackname}_Sum
  set is_movie = "no"
  ${proc_2dx}/linblock "Using aligned average ${input_image}.mrc as input file for CTF measurement."
# else
#   ${proc_2dx}/protest "ERROR: Aligned average ${movie_stackname}_Sum.mrc does not exist."
endif
#
if ( ${CS} == "ScriptWillPutNumberHere" ) then
  set CS = ${Default_CS}
  echo "set CS = ${CS}" >> LOGS/${scriptname}.results
endif
#
if ( ${KV} == "ScriptWillPutNumberHere" ) then
  set KV = ${Default_KV}
  echo "set KV = ${KV}" >> LOGS/${scriptname}.results
endif
#
if ( ${defocus_phase_shift_doit} == "-" ) then
  set defocus_phase_shift_doit = ${Default_phase_shift_doit}
  echo "set defocus_phase_shift_doit = ${defocus_phase_shift_doit}" >> LOGS/${scriptname}.results
endif
#
# if ( ! -e ${input_image}.mrcs ) then
#   ${proc_2dx}/protest "ERROR: Drift-corrected stack ${input_image}.mrcs does not exist."
# endif
#
if ( -e ${movie_stackname}_Sum.mrc ) then
  echo "# IMAGE-IMPORTANT: ${movie_stackname}.mrcs <DriftCor stack>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: ${movie_stackname}_Sum.mrc <DriftCor image (2D, no DW)>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: ${movie_stackname}_Sum_fft.mrc <DriftCor image FFT (2D, no DW)>" >> LOGS/${scriptname}.results
  echo "# IMAGE: ${movie_stackname}.mrc <DriftCor image (2D, with DW)>" >> LOGS/${scriptname}.results
  echo "# IMAGE: ${movie_stackname}_fft.mrc <DriftCor image FFT (2D, with DW)>" >> LOGS/${scriptname}.results
else
  echo "# IMAGE-IMPORTANT: ${movie_stackname}.mrc <DriftCor image (2D, with DW)>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: ${movie_stackname}_fft.mrc <DriftCor image FFT (2D, with DW)>" >> LOGS/${scriptname}.results
endif
#
##########################################################################
${proc_2dx}/linblock "Calling CTFFIND4.1"
##########################################################################
#
#set sample_pixel = `header ${input_image}.mrcs | grep spacing | cut -d' ' -f7`
if ( ${is_movie} == "yes" ) then
	set sample_pixel = `clip info ${input_image}.mrcs | grep Scale | cut -d'x' -f2 | tr -d " "`
else
	set sample_pixel = `clip info ${input_image}.mrc | grep Scale | cut -d'x' -f2 | tr -d " "`
endif

#
echo ":: "
echo "::Running:"
echo ":: "
echo "::"${app_ctffind} \<\<eot 
if ( ${is_movie} == "yes" ) then
	echo "::"${input_image}.mrcs
	echo "::"${is_movie}
	echo "::"4
else
	echo "::"${input_image}.mrc
endif
echo "::"diagnostic_output.mrc
echo "::"0.806 #${sample_pixel}
echo "::"${KV}
echo "::"${CS}
echo "::"${ampcon}
echo "::"1024
echo "::"${CTFFIND4_defocus_res_min}
echo "::"${CTFFIND4_defocus_res_max}
echo "::"${df_start}
echo "::"${df_end} 
echo "::"${df_step}
echo "::"no
echo "::"no
echo "::"yes
echo "::"1000.0
if ( ${defocus_phase_shift_doit} == "y" ) then
  echo "::"yes
  set CTFFIND4_defocus_phase_shift_L_rad = `echo ${CTFFIND4_defocus_phase_shift_L} | awk '{ s = $1 * 3.141592654 / 180.0 } END { print s }'`
  set CTFFIND4_defocus_phase_shift_H_rad = `echo ${CTFFIND4_defocus_phase_shift_H} | awk '{ s = $1 * 3.141592654 / 180.0 } END { print s }'`
  set CTFFIND4_defocus_phase_shift_step_rad = `echo ${CTFFIND4_defocus_phase_shift_step} | awk '{ s = $1 * 3.141592654 / 180.0 } END { print s }'`
  echo "::"${CTFFIND4_defocus_phase_shift_L_rad}
  echo "::"${CTFFIND4_defocus_phase_shift_H_rad}
  echo "::"${CTFFIND4_defocus_phase_shift_step_rad}
else
  echo "::"no
endif
echo "::"no
echo "::"eot
echo ":: "
#

echo "#IMAGE-IMPORTANT: diagnostic_output.txt <Summary of results (TXT)>" >> LOGS/${scriptname}.results
# echo "#IMAGE-IMPORTANT: diagnostic_output.mrc <Diagnostic images  (MRC)>" >> LOGS/${scriptname}.results
echo "#IMAGE-IMPORTANT: diagnostic_output_avrot.txt <Detailed results (TXT)>" >> LOGS/${scriptname}.results
#
echo "<<@progress: 20>>"
#
# set CTFFIND4_defocus_phase_shift_L_rad
# set CTFFIND4_defocus_phase_shift_H_rad
# set CTFFIND4_defocus_phase_shift_step_rad
if ( ${is_movie} == "yes" ) then 
if ( ${defocus_phase_shift_doit} == "y" ) then
time nice nohup ${app_ctffind}  << eot 
${input_image}.mrcs
${is_movie}
4
diagnostic_output.mrc
${sample_pixel}
${KV}
${CS}
${ampcon}
1024
${CTFFIND4_defocus_res_min}
${CTFFIND4_defocus_res_max}
${df_start}
${df_end} 
${df_step}
no
no
yes
1000.0
yes
${CTFFIND4_defocus_phase_shift_L_rad}
${CTFFIND4_defocus_phase_shift_H_rad}
${CTFFIND4_defocus_phase_shift_step_rad}
no
eot
#
else
#
time nice nohup ${app_ctffind}  << eot 
${input_image}.mrcs
${is_movie}
4
diagnostic_output.mrc
${sample_pixel}
${KV}
${CS}
${ampcon}
1024
${CTFFIND4_defocus_res_min}
${CTFFIND4_defocus_res_max}
${df_start}
${df_end} 
${df_step}
no
no
yes
1000.0
no
no
eot
#
endif
else
if ( ${defocus_phase_shift_doit} == "y" ) then
time nice nohup ${app_ctffind}  << eot 
${input_image}.mrc
diagnostic_output.mrc
${sample_pixel}
${KV}
${CS}
${ampcon}
1024
${CTFFIND4_defocus_res_min}
${CTFFIND4_defocus_res_max}
${df_start}
${df_end} 
${df_step}
no
no
yes
1000.0
yes
${CTFFIND4_defocus_phase_shift_L_rad}
${CTFFIND4_defocus_phase_shift_H_rad}
${CTFFIND4_defocus_phase_shift_step_rad}
no
eot
#
else
#
time nice nohup ${app_ctffind}  << eot 
${input_image}.mrc
diagnostic_output.mrc
${sample_pixel}
${KV}
${CS}
${ampcon}
1024
${CTFFIND4_defocus_res_min}
${CTFFIND4_defocus_res_max}
${df_start}
${df_end} 
${df_step}
no
no
yes
1000.0
no
no
eot
#
endif
endif
#
echo `tail -n 1 diagnostic_output.txt` | cut -d\  -f2-4 | sed 's/ /,/g' > tmp.1
set defocus = `cat tmp.1`
echo `tail -n 1 diagnostic_output.txt` | cut -d\  -f5 > tmp.1
set CTFFIND4_phase_shift = `cat tmp.1 | awk '{ s = $1 * 180.0 / 3.141592654 } END { print s }'`
set defocus_phase_shift = ${CTFFIND4_phase_shift}
echo `tail -n 1 diagnostic_output.txt` | cut -d\  -f6 > tmp.1
set CTFFIND4_CCvalue = `cat tmp.1`
set defocus_CCvalue = ${CTFFIND4_CCvalue}
echo `tail -n 1 diagnostic_output.txt` | cut -d\  -f7 > tmp.1
set CTFFIND4_RESMAX = `cat tmp.1`
set defocus_RESMAX = ${CTFFIND4_RESMAX}
\rm tmp.1
echo "set defocus = ${defocus}" >> LOGS/${scriptname}.results
echo "set CTFFIND4_phase_shift = ${CTFFIND4_phase_shift}" >> LOGS/${scriptname}.results
echo "set defocus_phase_shift = ${defocus_phase_shift}" >> LOGS/${scriptname}.results
echo "set CTFFIND4_CCvalue = ${CTFFIND4_CCvalue}" >> LOGS/${scriptname}.results
echo "set defocus_CCvalue = ${defocus_CCvalue}" >> LOGS/${scriptname}.results
echo "set CTFFIND4_RESMAX = ${CTFFIND4_RESMAX}" >> LOGS/${scriptname}.results
echo "set defocus_RESMAX = ${defocus_RESMAX}" >> LOGS/${scriptname}.results

#
echo "::Estimated resolution limit by EPA: ${CTFFIND4_RESMAX}"
#
set tmpdef1 = `echo ${defocus} | sed 's/,/ /g' | awk '{ s = $1 } END { print s }'`
set tmpdef2 = `echo ${defocus} | sed 's/,/ /g' | awk '{ s = $2 } END { print s }'`
set defocus_angle = `echo ${defocus} | sed 's/,/ /g' | awk '{ s = $3 } END { print s }'`
set CTFFIND4_defocus = `echo "scale=3; ( ${tmpdef1} + ${tmpdef2} ) / 20000.0 " | bc `
set defocus_defocus = ${CTFFIND4_defocus}
set tmp = `echo "scale=3; sqrt((( ${tmpdef1} - ${tmpdef2} ) / 2)^2) / 10000.0 " | bc `
set defocus_astig =  `echo ${tmp} | awk '{ s = $1 + 0.0 } END { print s }'`
set CTFFIND4_defocus = ${defocus_defocus}
echo "::Average defocus = ${CTFFIND4_defocus} microns"
echo "set CTFFIND4_defocus = ${CTFFIND4_defocus}" >> LOGS/${scriptname}.results
echo "set defocus_defocus = ${defocus_defocus}" >> LOGS/${scriptname}.results
echo "set defocus_angle = ${defocus_angle}" >> LOGS/${scriptname}.results
echo "set defocus_astig = ${defocus_astig}" >> LOGS/${scriptname}.results
#
#################################################################################
${proc_2dx}/linblock "Running: labelh.exe to normalize MRC diagnostic image"
#################################################################################
# 
echo "<<@progress: 90>>"
#
\rm -f CTFDiag.mrc
#
time ${bin_2dx}/labelh.exe << eot
diagnostic_output.mrc
42
CTFDiag.mrc
eot
#
echo "#IMAGE-IMPORTANT: CTFDiag.mrc <CTF Diagnostic Plot (MRC)>" >> LOGS/${scriptname}.results
#
##########################################################################
${proc_2dx}/linblock "Update statuspage images."
##########################################################################
#
\rm -f tmp.png
\rm -f tmp2.png
\rm -f tmp3.png
\rm -f CTFDiag.mrc.png
\rm -f STATUS/4-image.jpg
${app_2dx_mrc_converter} --size 400 CTFDiag.mrc tmp.png
${app_python} ${proc_2dx}/PNGannotator.py tmp.png tmp2.png 10 345 0 "CTFFIND4.1 Thon ring fit"
${app_python} ${proc_2dx}/PNGannotator.py tmp2.png tmp3.png 10 360 0 "${date_text}"
${app_python} ${proc_2dx}/PNGannotator.py tmp3.png CTFDiag.mrc.png 10 375 0 "Defocus: ${CTFFIND4_defocus} um.  CTF Resolution: ${CTFFIND4_RESMAX} A"
${app_python} ${proc_2dx}/PNGannotator.py tmp3.png STATUS/4-image.jpg 10 375 0 "Defocus: ${CTFFIND4_defocus} um.  CTF Resolution: ${CTFFIND4_RESMAX} A"
\rm -f tmp.png
\rm -f tmp2.png
\rm -f tmp3.png
#
echo "<<@progress: 100>>"
echo "<<@evaluate>>"
#
##########################################################################
${proc_2dx}/linblock "${scriptname} - normal end."
##########################################################################
#
exit
#
